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We present a new stellar evolution code and a set of results, demonstrating 
its capability at calculating full evolutionary tracks for a wide range of masses 
\ and metallicities. The code is fast and efficient, and is capable of following 

t*^* , through all evolutionary phases, without interruption or human intervention. 

fNj ■ It is meant to be used also in the context of modeling the evolution of dense 
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stellar systems, for performing live calculations for both normal star models 



ON i 

' and merger-products 

OC ' 

C^) . The code is based on a fully implicit, adaptive-grid numerical scheme 

that solves simultaneously for structure, mesh and chemical composition. Full 
' details are given for the treatment of convection, equation of state, opacity, 

h ; 

nuclear reactions and mass loss. 

Results of evolutionary calculations are shown for a solar model that 
matches the characteristics of the present sun to an accuracy of better than 
1%; a 1 M model for a wide range of metallicities; a series of models of 
stellar populations I and II, for the mass range 0.25 to 64M©, followed from 
pre-main-sequence to a cool white dwarf or core collapse. An initial final- 
mass relationship is derived and compared with previous studies. Finally, we 
briefly address the evolution of non-canonical configurations, merger-products 
of low-mass main-sequence parents. 
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2 A. Kovetz, O. Yaron and D. Prialnik 
1 INTRODUCTION 



Simulating the evolution of a star requires the solution of a set of partial differential equations 
with boundary conditions at the center and surface, involving extensive input physics, such 
as equations of state, nuclear reactions, opacities, as well as recipes for treating convection, 
mass loss, or material mixing. This is accomplished by complex computer codes that are 
time consuming and depend on a large number of adjustable parameters, both physical and 
numerical, needed for dealing with evolutionary phases that are different in nature. From 
the formation to the death of a star, differences between evolutionary phases are so large, 
that studies are usually devoted to — and often codes are devised for — a specific part of a 
star's life, ignoring or simplifying, or suppressing others. So far, no code has been suited or 
applied to obtain complete, unabridged evolutionary tracks over the entire range of stellar 



masses and metallicities, a 



Pols et al. 



1995, 



Pols et al- 



thoug h many have come close to accomplishing this task (e.g. 



19981 ). For example, most (if not all) evolution codes crash at 



the helium core flash phase. Most of the stellar evolution codes do not solve simultaneously 
for the structure and the composition; this introduces serious errors in some critical phases 
whenever the mass grid {rrii} changes, as it must eventually (jStancliffd 120061 ) . Our aim has 
been to develop a versatile and robust stellar evolution code that is free of such handicaps. 

A further demand on the code is efficiency and speed. Furthermore, it should be capable 
not only of evolving any star through all phases without intervention, but also of dealing 
with peculiar objects. Such a code could be incorporated into an N-body code that deals 
with dense stellar systems, if not at present, then — given the rapid and continual advance 
in computing power — in the foreseeable future. The computation methods of N-bo dy grav- 



i tating systems have u ndergone a rev olutionary development owing to the work of 



Aarseth 



( 1963 ) (see review 



Heggie fc Hut 



2003 



Aarseth 



1999) a nd gaining impetu s in the past two decades (e.g. 



Hurley et al. 



2001 



Hurley et al. 



20051 ): not onl y have new algorithms 



been developed, capable of dealin g with dense stellar systems (e.g. iPortegies Zwart et al 



2001 



Portegies Zwart et al.ll2004l ). bu t also special 



rardware has been constructed under 



the GRAPE (GRAvity PipE) project flMakino et al 



19971 ) 



However, in order to render these sophisticated N-body calculations realistic, the effect 
of the structure and evolution of the constituent stars must be considered as well. This led, 
less than a decade ago, to the development of the MODEST (Modelling DEnse STellar 
systems) project, whose aim is to combine N-body dynamics with the hydrodynamics of 
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stellar collisio ns on the one han d, and with stellar evolution of the cluster population, on 



the other (see 



HutetaL 



(120031 )). So far, studies of stellar systems have resorted to short 
cuts based on sets of discrete pre-calculated evolutionary tracks: either interpolating between 
them, or using parametrized fit formulae. Clearly, this procedure is incapable of dealing with 
'non-canonical' stars, the outcome of collisions and mergers. 

In this paper we thus present a new evolutionary code that we have developed with this 
aim in mind. The outline of the code and method of solution are presented in the next 
section, Section EJ the input physics is described in some detail in Section [3l and results of 
representative calculations are discussed in Section HI 



2 THE EVOLUTION CODE 

2.1 Set of Equations and Boundary Conditions 

The equations that govern the evolution of a star are those of continuity, hydrostatic equi- 
librium, energy transfer (radiative or convective), energy balance, and composition balance: 



d 4vr 

dm 3 

dp 

dm 
d\nT 

dm 
du 

di +P dip~ 



P 

Gm 
4 7rr 4 



7 <91np 

dm 



d 1 
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dm 



BY, 



dYj 

3 dm 

dF 
R* - 3 



(1) 
(2) 
(3) 
(4) 
(5) 

dt ^ J dm ^ 
In these equations, mass m and time t are the independent variables. The dependent 

ones are radius r, density p, temperature T, and the number fractions Yj, related to the 

mass fractions Xj by Yj = Xj/Aj, where Aj is the j'th atomic mass. The particle flux Fj of 

the j'th species is assumed to be diffusive (proportional to the abundance gradient of the 

j'th species), determined by the diffusion coefficient aj. 

We regard (p, T, Y) as the basic thermodynamic variables. They determine, through the 

equation of state, the pressure p(p,T,Y) and the specific energy u(p,T,Y), as well as the 

opacity k(p, T, Y) , energy production rate q(p, T, Y) and nuclear energy rates Rj (p, T, Y) (via 

an imported list of tables and formulae). The temperature 'gradient' V(r, L, m, p, T, Y) and 
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the (convective mixing) diffusion coefficients <jj(r, L,m, p,T,Y) are provided by convection 
recipes. 

The foregoing equations are to be solved subject to the following boundary conditions: 
at the centre, 



r = 0, L = 0, Fj = 
at the surface, which we take to be the stellar photosphere, 



Kp G 



[l-T)gT s , L = Anr 2 aT 4 



Fj = 



(7) 



(8) 



In the first member of eq. (jHJ), pc is the material ('gas') pressure which, together with the 
radiation pressure pn, makes up the total p\ Y = KL /AircGm,-, q — Gm/r 2 : and r s is the 



photospheric optical depth, which we take to be unity (IKovetz 



1998 . 



Kovetz 



1999|). 



We shall solve the e quations o 



M, but we shall follow (lEggleton 



evolution over a grid of mass points mi = 0, rri2, • • • , m n 



1971 



Eggletonlll972l ) in using an adaptive grid, where the 



mass points m 2 , . . . ,m„-i depend on the solution. Since, by eqs. ([I])-©, r at the centre 
varies like m 1//3 , and p like m 2 / 3 , we replace m by x = m 2//3 , and r by s = r 2 , in this pair of 
equations. Equations (CD)-© then become 

3 rx\ I 

3G /£\2 

s. 



ds 
dlnp 



QTCp 

dlnT = Vdlnp 



© 2 - 



dL 



q 



5u + pS- 



5t 



dm 



-a 



dm 

SY : 

Rj r^- 1 dm 

ot 



(9a) 
(9b) 
(9c) 
(9d) 

(9e) 
(9f) 



These may be regarded as differential equations, written in terms of differentials; alterna- 
tively, they may be thought of as representing difference equations. In the latter case, at 
the centre, the indeterminate ratio x/s = 0/0 is replaced by its limit (47rp 1 )/3) 2/ ' 3 , where p\ 
is the central value of the density. The change from r to s obviously requires appropriate 
changes (such as s = 0, g = Gm/s) in the boundary conditions. 

The equations of structure and composition are solved simultaneously with a mass dis- 
tribution function, implementing an adaptive mesh. This is done by requiring constant in- 
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crements of a monotonic function of the form 



/ = (m/M) 2/3 + dX H - c 2 lnp - c 3 In • 



T 



tw ( 10 » 

where the c's are appropriate non-negative constants. Near the centre the requirement of 
equal increments of / will lead to equal increments of x = m 2 / 3 . The second term of / will 
force equal increments of the hydrogen mass fraction where Xjj changes rapidly (at an H- 
burning shell). The third will lead to equal steps of hip towards the surface, where m/M w 1 
and Xh is uniform; and the last term will cause a fine subdivision around T = C4 ~ 20, OOOfT, 
where the opacity varies rapidly over several orders of magnitude. 



2.2 Numerical Scheme 

The variables (s, m, L, p, T, p, Yj) are represented by arrays over a grid of % = 1, . . . , n, where 
i — 1 corresponds to the centre, and i — n to the surface (photosphere). Thus eqs. (J9al) -(l9bl) 
become the difference equations 



hxpi — hip 



1 
2 

i-l 



1/2 



3G 



+ 



2 



3G 



1/2" 



(Xj Xj_i) 



(11) 



;i2) 



There is one such pair of equations for each % = 2, . . . , n. Together with the boundary 
conditions s± = and {hpg)u = (1 — r n )g„r s , these add up to 2n equations. 

The variables V, L and Fj, related to the energy and particle fluxes, are replaced by 
arrays that refer to the midpoints i ±1/2. Thus eqs. (!9c|) - (l9f|) become 

InT; - lnTj_i = Vj_i(lnpj - lnp^i) 



- (T i+ 1 - 



8t 



Pi 



5- 

pi 

St 



Y{+1 — Yi 



2 rn i+ i 



mi 



(13) 
(14) 

(15) 
(16) 



where, in the last pair of equations, we have suppressed the index j that refers to the nuclear 
species. The coefficients Vj_i and a i+ i are evaluated by using the arithmetic means of the 
grid-point arguments, for example r € _i = (rj_i + ri)/2. Again, there is one eq. f)13p for 
each % = 2, . . . , n, which, together with the boundary condition L n = Airr^aT^, brings the 
number of equations up to 3n. Furthermore, there is one set of eqs. (f!4"|) and fflBT) for each 
% = 1, . . . , n (and for each one of the species). If J is the number of species, the number of 



6 A. Kovetz, O. Yaron and D. Prialnik 

equations becomes (4 + J)n. At i = 1 we set L { _i = F { _i = and rn,j_i = in eqs. ([T] 
and ( |T6i) . which takes care of the central boundary conditions L = F = m = 0. At i = n we 
set = L n , F i+ i = and m^x = m n in eqs. (JHJ) and (fl6l) . This is in accord with the 
surface boundary conditions. 

The requirement of equal increments of the mesh function / is simply 

fi+l — ft — fi ~ fi-l- (17) 

There is one such equation for each % = 2, ... ,n — 1. At the ends % = 1 and i = n we 
respectively impose the two boundary conditions 

mi = 0, m n = M + MSt, (18) 

where M(m n , r n , L n ) is the rate of mass accretion — or loss, if negative. Thus we have a total 
of (5 + J)n equations for (5 + J)n variables — 5 arrays (s, L, m, p, T) and J arrays Yj, each 
array being of length n. 

The partial time derivatives, du(m,t)/dt, etc., have been replaced, respectively, by dif- 
ference ratios 5u/5t, etc. When a configuration at a previous time is available, Su is usually 
taken to be u(m,t) — u(m,t — St), where u(m,t) is iterated upon. The solution of eqs. 
(J9a|) - (j9fj) then has the accuracy 0(6t). It should be noted that u(t — St) is represented by 
a grid function over a (previous) set of m^'s that will not generally include the m for which 
u(m,t — St) is desired. We therefore determine u(m,t — St) by interpolation, using cubic 
Hermite splines. These splines have the advantage that, if the grid function vanishes at two 
consecutive m^s, the interpolant will not dip below zero anywhere between them. This is 
especially important when interpolating the number fractions Yj. 

Except at the first time step, the previous, as well as the anteprevious, configurations are 
available. Instead of a chord through u(m,t) and u(m,t — St), we can then pass a parabola 
through u{m,t), u(m,t — St) and u(m,t — St — St'), and evaluate its derivative at t. If this 
derivative is again denoted by Su/St, we have 

Su = au(m, t) + f3u{m, t — St) + 7«(m, t — St — St'), (19) 

where 

_ St' + 2St R _J t ' + 5t _ ( st ) 2 

a ~ St' + St 1 P ~ 8V~> 7 ~ (St' + St)Sf ( } 

This leads to a solution with accuracy 0(St 2 ). Of course u(m, t — St — St'), like u(rn, t — St), 
has to be determined by interpolation. 

The (5 + J)n nonlinear eqs. ({IT]) - ({TBI for the the arrays (s, L,m, p,T,Yj) are solved 
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simultaneously by Newton-Raphson iterations. This requires, at each iteration stage, the 
solution of a linear system with a band matrix of order (5 + J)n, and band width 15 + 4J. 

The derivatives required by the Newton-Raphson method are evaluated analytically 
whenever possible. In the case of opacities, which are obtained from tables with the aid of 
cubic Hermite spline interpolation, we use the (analytic) derivatives of the splines. Numeri- 
cal derivatives are only used for the energy generation and loss rates, because the neutrino 
loss rates are provided by cumbersome fit formulae. 

2.3 Computational Details 

Our automatically varying timesteps, determined mainly by limits imposed on the maximal 
changes (a few percent), and on the number of Newton-Raphson iterations, allowed during a 
timestep, span a wide dynamic range — from seconds/minutes during core or shell flashes to 
several times 10 8 or even 10 9 years in the main-sequence phase (of low- mass stars). With a 
relative accuracy of ~0.0001, the typical number of Newton-Raphson iterations is 3-4. The 
grid mass shells, determined by the mass-distribution function, span a range of ~ 1O~ 15 M 
(in a WD atmosphere) to > 1CT 1 M (in an inert stellar core). There is an option of fixing 
the mass grid, which we are forced to use during the WD cooling phase, when the mass 
array {m^} ceases to be monotonically increasing in double precision arithmetic. With these 
features in mind, the typical number of grid points may be as low as 150 or 200; a typical 
number of timesteps for a complete evolutionary track is 1000; and typical execution time 
is of the order of 10 (±5) minutes on a portable computer (Pentium 4 and higher). The 
latter is, however, strongly dependent on both physical behaviour (e.g. mass-loss rate or 
the amount of evolutionary phases taking place) and computational prescriptions (required 
outputs / interfaces) . 

The code — targeted for Unix/Linux machines — is written in Fortran 90 and consists of 
an online graphical interface using Tim Pearson's PGPLOT. 

3 INPUT PHYSICS 

3.1 Equation of State (EOS) 

The EOS is derived from a free energy, which is a sum of ionic, radiative and electronic 
contributions, together with corrections for pressure ionization, Coulomb interactions and 
quantum effects: 
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F=j2 F ( T > v > N ^-l aTiv 

i 

+ fi(T,V> F ) + (up - mc 2 )N 
+ Fpi + Fcq. 



(21) 



where V is the volume, a is Stefan-Boltzmann's constant, fip is the Fermi chemical potential, 
and m is the electron's mass. The free energy of the iVj particles of the i'th ionic species is 



TV- 

F(T, V, Ni) = (fcTln — — kT + Xi )N u 

Zi 



Zi 



VQ 



3 



h 



\J1ixmikT ' 



(22) 



where Xi is the reference energy (relative to the completely ionized state) of the z'th ion, 
and Qi is its partition function. The thermal length l{ depends on the temperature and on 
the z'th particle's mass mj = AitriH, where A$ is the atomic or molecular weight. 

We take account of ionization equilibria for hydrogen and helium; heavier elements (the 
'metals') are assumed to be completely ionized. In the stellar envelope, where the metals 
amount to at most a few percent by mass, and a few thousandths by number, this introduces 
an error that is much smaller than other uncertainties in the EOS. In a carbon/oxygen stellar 
core, the metals are pressure-ionized in any case. Ionization equilibria of the metals play 
an important role in determining the opacity, but we use opacity tables that are entirely 
independent of our EOS. 

Remembering that the reference energies for the completely ionized species H + and He ++ 
are zero by definition, the Xi f° r H, H 2 , H + , He, He + and He ++ are, respectively, -13.598, 
-31.673, 0, -79.003, -54.416 and (in eV). Also, Xi = for the metals. 

Except for the case of H 2 , we replace the partition function Qi by a constant statistical 
weight gi, which is 1 for H + , He, He ++ and all metals, and 2 for H and He + . For the hydrogen 
molecule, w e use our own table of Qho(T ) , whi ch we have calculated, using the molecular 



constants of 



Tatuml (119661 ) ; see also llrwinl (119871 ) . 



Electrons and positrons are described by the fermion grand thermodynamic potential 



fi(T, V, pLp) = -Cmc 2 V I T(e/P)D+(e, 0) de/P, 



(23) 
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where 

c-^) 3 r ( x) = |(«-- 1,1, 

D±(e, d>) = — t ± — 7 . 

"=*r- * = < 24 > 

(IRakavy et al .11196 71 ) . The Fermi chemical potential /xp, which includes the rest-mass energy, 
is connected with the number difference, electrons minus positrons, through 

N = -Q flF (T,V,fi F ) 

POO 

= CV / r 7 (e/^)D_(e, 0) de//3 = iV e - N p , (25) 

J/3 

where Q^ F (T, V, fip) denotes the partial derivative dfl(T,V, /xp) jd\ip. Clearly the positron 
contribution, which is due to the second term of D±, becomes insignificant whenever cb is 
large (say <fi ^ 15). The last equation determines Hf(T, V, N) as a function of T, V and N 
(actually the Fermi parameter ep = (fj,p — mc 2 )/kT in terms of T and N/V). The number 
difference N must satisfy the equation of charge neutrality 

N = Y,Z l N l . (26) 

If, in the expression (1211) for the free energy, electrons appeared only in the second line, 
then it would follow that 



2 



F N (T, V, N) = fi MF /iF,Af + ^f,nN + (jlf- mc 

= jj,p — mc 2 , (27) 

where the subscript iV denotes the partial derivative with respect to N, at constant T and 
V. Thus up would indeed be the electron chemical potential \i = F N + mc 2 . We maintain 
the distinction (between /i and fip) because other parts of the free energy — for example the 
pressure ionization term Fpj — too, depend on the electron number density. 

The pressure p = —Fy, the entropy S = —Fp and their derivatives require derivatives of 
Q, with respect to (3 or <f), up to the second order. This leads to five additional Fermi-Dirac 
integrals, in which T(e//3) is replaced by T'(e/(3), (e/ /3)T'(e/ '/?), T"(e/p), (e/ (3)T"(e/ (3) or 
(e//5) 2 r"(e//5). In the degenerate case, when ep = (ft — (3 > 5, Q is calculated by Sommerfeld's 
method, and then differentiated. Otherwise the six integrals involving the first, electronic, 
part of D± are calculated in one swoop, using Gaussian quadrature. The nodes and weights 
for this quadrature are calculated at the beginning of the run, and their number can be 
chosen by the user (the code's default is 12 nodes). The positronic contribution, which is 
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due to the second part of D±, is then obtained by using the same procedure, with <fi replaced 
by —0. The last step is only carried out when < 15; otherwise, positrons a re ign ored. 



Pols et al. 



(|1995h : 



" (ci/a ' )C2 b + eF + c 3 ln(l + x/c 4 )], 
N e /V, y = 13.60/fcr, 



The pressure ionization term in the free energy is taken from 

Fpi = -N e kTg(n e , T) + N e0 kTg(n e0 , T), 
g(n e ,T) = e 
x = n e rriH, n e 

(ci,c 2 ,c 3 ,c 4 ) = (3,0.25,2,0.03), (28) 

where the 13.60 is in eV, and the units of C\ and c 4 are g cm -3 . Furthermore, N e q is the total 
number of electrons, bound or free, and n e0 = N e0 /V. The object of Fpi is to induce pressure 
ionization by reducing the electronic chemical potential as the number of electrons n e a^ in 
a cube with side do, the Bohr radius, increases. Of course Fpi tends to zero as ionization 
becomes complete, that is, as N e — > N e0 . 

The last term, Fqq, in the free energy depends on the Coulomb parameter T and on the 
Debye parameter A. For a one-component plasma (OCP) 

Z?e 2 



TjkT ' 



where Zj is the atomic number, e is the electron charge and = 
ion-sphere radius. For a mixture, we replace this by 

^XiZf/A \(Y, x iZi/ A i\ 2 ^ n eY ,z f_ 

kT' 



(29) 

(A-nNi/W)- 1 ^ is the 
(30) 



Y.XiZi/Ai LV Y.Xi/Ai ) 3 
where the sums refer to a fully ionized plasma mixture with mass fractions X^. Again, for a 

one component plasma, 



A; 



pi 



4irZ?e 2 ni 



kT 



to, 



pi 



rrii 



where ui pi is the plasma frequency. For a mixture, we replace this by 



A 



fUjJr, 



n,, 



(31) 



(32) 



kT p 

where A^4 is Avogadro's number. The expression for Fqcp t akes different form s for the gas- 



lb en et al 



(|1992f ). Noting that 



liquid and for the solid phases, and is based on the work of 
the OCP form of the translational part (that is, setting Qi = 1 and omitting the Xi' s ) °f the 
ionic free energy ^ F(T, V, Nj) is 



F° 1 7T 

— Ar ,^ = 3 In A - 1.51nr + - In - - 1 

£ NikT 2 6 

= 31nA - 1.51nr - 1.32351, 



(33) 
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Iben et al. 



(1992) write the OCP free energies in the form 
Fqcp 



Fqcp 



3 In A - 1.5 In T - 1.32351 - H^T) + J(A), 

where J (A) takes care of the quantum effects in the gas-liquid phase, 

Hi =— r 3/2 + T 3 (-0. 104584 
3 

+ 0.172110 lnT -0.033724r 3/2 ), T ^ 1, 

H t =0.897744r - 3.801720r 1/4 + 0.758240r~ 1/4 

+ 0.814871 lnT + 2.584778, 1 ^ T 200, 

1612.5 



(34) 
(35) 



H s =0.895929r + 



r 2 



(36) 



and .Fvih(A ) 



Iben et al. 



is t he vibrational contribution to the free energy (IKovetz fc Shavivl Il970l ). 



(Il992l ) have shown that F vi b/ Yl AfikT can be fitted by a weighted sum of two 



Debye free energies: 



vib 



t A s , i A \ 

x _ =aL(— ) + (l-a)L(— ). 



where a = 0.5711, Ai = 1.0643, A 2 = 2.9438, and L(x) is given by 



L{x) = lx + 3 ln(l - e _x ) - D(x) 



x c 



t 3 dt 
e* - 1" 



(37) 



(38) 



The function J (A) is known (IShaviv fc Kovetzlll972l ) to have the the high-temperature limit 
A 2 /12. At low temperatures the OCP liquid should resemble a bec lattice, with the ions 
vi brating abo u t thei r equilibrium positions. This leads to a J (A) proportional to A: according 



to 



Iben et al 



(11992( 1. J (A) -> 1.06980 A (although their foregoing fit for F vib yields F vib -> 
0.76758A). They then suggest a functional form for J (A) that interpolates between these 
limits. But this leads to a non-monotonic entropy (T-derivative of the gas-liquid -Focp); i n 
particular, the specific heat has the required T 3 dependence at low T, but with the wrong 
sign! 

Rather than adopt Iben et al.'s J (A), we note that, for A << 1, F vib /^A^fcT tends to 
3 In A — 1 — 1.49602, whereas for A >> 1 it tends to 0.76758A, and therefore set 



OCP 



vib 



+ 1.49602 - 0.32351 - 1.5 In T - Hi(T). 



(39) 



^NikTJug ZNikT 
The OCP free energies include the contribution of the translational degrees of freedom. Since 
our free energy already includes ^2 F(T,V, Ni), we must, in order to obtain Fqq, subtract 
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F° from each one of the Focp's. Thus, finally, 

3 lnA + 2.49602 -ifi(r) 



FcQ \ -^vib 



ZNikTJn q J2NikT 

(ewL = A " 31nA + L32351 + L5 lnr - ^( r )- ^ 

Formally, the difference between the liquid and solid free energies leads to a phase transition 
when 

1.49602 -^(r) = 0.32351 + 1.5 In T - H 8 (T). (41) 

The root of this equation, the 'melting T\ is r m = 178.2119 . We avoid this complication by 
interpolating for Fqq in the interval (T m — 2, T m + 2). 

We shall not pause to write down the equations — such as pn 2 — for H 2 <-> 2H, or 
^He = ^He+ + U ~ rac 2 for H e <-» He + +e — that determine the various states of hydrogen or 



helium (e.g. iPols et all (11995h ). 



3.2 Opacities 

The opacities, which generally depend on density, temperature and composition, are of two 
kinds: radiative and conductive. For the radiative part we use Boothro yd's interpolation 



progr a to interpolate within the OPAL Rosseland mean opacity tables ( llglesias fc Rogers! 



19961 ). Each one of the OPAL tables is for a given hydrogen mass fraction X, a given total 



heavy element mass fraction Z (distributed in accordance with one of a number of standard 
'mixes'), a given carbon mass fraction excess Xq (such that the total carbon mass fraction 
is Xc , plus the carbon mass fraction contained in Z) , and a given oxygen mass fraction 
excess Xo ■ The helium mass fraction is of course 1 — X — Z — Xc — Xo ■ 

Each one of the OPAL tables spans a temperature range 3.75 < logT < 8.70 and a 
range —8 < logi? < +1 of logi? values, where R = p/T§ , with a cutout at the high T, high 
R corner, and sometimes at the low T, low R, corner. Boothroyd's interpolation program 
provides the OPAL opacity k, together with its density and temperature derivatives. (In 
this section, T is in degrees Kelvin, p in gr cm -3 , and k in cm 2 g _1 .) 



At the low temperature end the OPAL opacities are supplemented by the lFerguson et al 



( 120051 ) tables. These span a temperature range 2.70 < logT < 4.50, and the same R range 



as the OPAL tables. But their Z range has the upper limit Z = 0.10, and there is no 



A Stellar Evolution Code for Calculating Complete Tracks 13 
provision for C or O excesses. We interpolate among them with a Z value equal to the lesser 
of Z + Xc + Xo and 0.1. 

At the high temperature end, logT > 8.70, we exte nd the OPA L opacities by using 



electron/positron scattering opacity according to the fit of llbenl (119751 ): 

K es = [0.2 -D-(D 2 + 0.0004) 1/2 ]2n ep /(iV A p), D = 0.05(logT 6 - 1.7), (42) 
where n ep is the sum of the electron and positr on number densities. 



Cassisi et al. 



(120071 ) tables. These span the 



Electronic conductivities are taken from the 
temperature range 3 < logT < 9 K, and the density range —6 < logp < 9.75 . There is one 
such table for each value of the atomic number Z ion , in fact 15 tables spanning the range 



1 < Z ion < 60. We use the interpolation program provided by lCassisi et al.l ( 120071 ). with Z ion 
equal to the square root of the average (by number) squared atomic number 

(^zfx./A.y^Xi/A). 

The conductivity is converted to a conductive opacity and — harmonically — combined with 
the radiative opacity. 

The various opacity interpolation programs provide the opacity n, together with its 
density and temperature derivatives. But an evolution code that simultaneously solves for 
the stellar structure and composition requires the derivatives of k with respect to composition 
as well. One way to get these is to evaluate the opacity at neighbouring compositions and 
then form difference ratios. 

Alternatively, we use the following method, which yields continuous opacity derivatives: 
at the beginning of the evolutionary run, we use the various interpolation programs to create 
a set of total — radiative and conductive — opacity tables that, for the initial stellar model's 
Z, span the triangular region of Fig. [TJ Along the x-axis of this figure we have seven values 
of the hydrogen mass fraction X, from to 1 — Z , with no carbon or oxygen mass excesses. 
Along the y-axis there are seven values of the combined C/O excess X C o — X c + X , 
again from to 1 — Z . Each point with positive Xco corresponds to a pair of tables: one 
with carbon excess equal Xco and zero oxygen excess (that is, 'excess all carbon'), and the 
other one with the same total excess Xco, but 'excess all oxygen'. The Xco > tables 
have the lower limit logT = 4.00, because the low-temperature Ferguson- Alexander tables 
correspond to zero C/O excesses. 

There are no points to the right of the hypotenuse (because they would correspond to 
negative helium mass fraction 1 — Z — X — Xco )• Our total number of opacity tables is 49, 
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and these replace the much larger number of OPAL, Ferguson- Alexander, and Cassisi tables. 
During MS hydrogen burning, the stellar core follows a path — from right to left — along the 
x-axis. During the HB (Horizontal Branch, core helium burning), C/O excesses rise and the 
core follows an upward path along the y-axis. Convective mixing may require the evaluation 
of the opacity in material containing both hydrogen and C/O excesses, that is, at points 
inside the triangle of Fig. [TJ 

In interpolating within the set of tables represented in Fig. (TJ we distinguish between 
three cases: 

Case I — No CO excesses. Interpolation is performed only within the 7 x-axis (hydrogen) 
tables. Within each table we use cubic Hermite splines to interpolate in logi? and logT, 
in order to obtain \ogn and its T and R derivatives. Among the seven resulting values of 
logK, we then interpolate in order to obtain the final opacity value, together with its X 
derivative, for the required hydrogen mass fraction. Similar interpolations among the seven 
T derivatives, and among the seven R derivatives, yield the T and R derivatives for the 
required X . (Since R— p/T 6 3 , the p derivative is simply related to the T and R derivatives.) 

Case II — no hydrogen — interpolation within y-axis tables (where for each Xco value 
there are two tables, the excess being completely in C for one, and completely in O for the 
second). We begin as in case I, by interpolating first among the 'excess all carbon' tables, 
and then among the 'excess all oxygen' tables. The final value of log k is then obtained by 
linear interpolation: 

log/t = ^-logK C + -^-hgno, (43) 
^co -^co 

logKc*, logKo denoting the opacities as obtained separately from tables for which excesses 
are all in C and and from tables for which excesses are all in O, respectively. Composition 
derivatives of the opacity, with respect to Xc, or with respect to Xq, are then obtained from 
the last formula. The final T and R derivatives are obtained by similar, linear interpolations. 

In comparing log ft obtained by this method with the one returned by Boothroyd's in- 
terpolation (which has its own uncertainties), we found deviations of no more than a few 
percent. And the largest of these were at fairly low temperatures, ~ 5.5 < logT < 6.5, 
where CO-rich opacities are less likely to be needed. 

Case III — both C/O excess and hydrogen — interpolations inside the triangle of Fig. [H 
This is a combination of Cases I and II. 

In Fig. [2] we display opacity profiles, alongside temperature, density and composition 
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profiles, at three snapshots during a solar model evolution ( §4.11) — Mid-MS, tip of RGB and 
the end-state as a cooling WD. Note that ranges of the y-axes values differ in between the 
three snapshots (columns), and it is apparent that there is a general decrease in opacity 
with the advance of evolution. The peak of opacity at low temperatures (around a few 10 4 
K), close to the surface, is due to the ionization of hydrogen. 



3.3 Nuclear Reaction Rates 



We use the following nuclear reaction network: 



1 H(p,/3+i/) 2 D(p, 7 ) 3 He 

3 He( 3 He,2p) 4 He 

3 He( 4 He,7) 7 Be 

7 Be(e-,z/) 7 Li(p,a) 4 He 

7 Be(p, 7 ) 8 B(/?+z/) 8 Be*(a) 4 He 

12 C(p, 7 ) 13 N(/3 + z/) 13 C(p, 7 ) 14 N 
14 N(p, 7 ) 15 0(/3 + z/) 15 N(p,a) 12 C 
14 N(p,7) 15 0(/? + z/) 15 N(p, 7 ) 16 
16 0(p, 7 ) 17 F(/5+z/) 17 0(p,a) 14 N 



4 He(a) 8 Be*(a, 7 ) 12 C 

12 C(«, 7 ) 16 

14 N(a,7) 18 F(|a,7) 20 Ne 

16 O(a, 7 ) 20 Ne 

20 Ne(a,7) 24 Mg 

12 C( 12 C,a) 20 Ne 

16 0( 16 0, 7 ) 32 S(7,a) 28 Si 

24 Mg(a, 7 ) 28 Si 

20 Ne( 7 ,a) 16 O 

24 Mg( 7 , a) 20 Ne 



with rates taken from 



The enhancement of the nuclear reactions 



Graboske et al 



Caughlan fc Fowlerl (1198 
by ele ctron screening is taken into account by following the prescriptions of 
(119731 ). 

Where several reactions are written in a chain, the later reactions are taken to be in 
transient equilibrium with the first one. The first five reactions — which constitute the pp- 
chain — are also assumed to be in transient equilibrium with each other, so that only the 
two major isotopes, 1 H and 4 He, need to be followed. Similarly, in the next four reaction 
chains — which constitute the CNO cycle — only the major isotopes 12 C, 14 N, and 16 are 
followed, and all other isotopes are taken to be in transient equilibrium. 

The triple-alpha reaction 4 He(a) 8 Be*(a, 7) 12 C, together with the four following lines, 
constitute helium burning, which involves two further major isotopes— 20 Ne and 24 Mg. The 



2 Website http://www.phy.ornl.gov/astrophysics/data/cf88/ . 



16 A. Kovetz, O. Yaron and D. Prialnik 



reaction 18 F( ^a, 7) 20 Ne is of course a fiction (IPols et al.lll995l ). intended to avoid the creation 
of 22 Ne, which is thus replaced by 20 Ne. 

Carbon burning proceeds — with comparable probabilities — through the two main branches 
12 C( 12 C,p) 23 Na , 12 C( 12 C,a) 20 Ne. Since the protons released by the first one interact with 
other species, in particular through the reaction 23 Na(p,a) 20 Ne , the net re sult of carbon 
burning can be described by the single reaction 12 C( 12 C,a) 20 Ne (jlliadisi 120071 ). 

Oxygen burning proceeds via many branches: the main product is 28 Si, with 32 S a close 
second {ibid.). We take 28 Si as our last major isotope. Thus, after oxygen burning, our 28 Si 
mass fraction is actually the sum of X( 28 Si) and X( 32 S). 

Carbon burning 12 C( 12 C,a) 20 Ne, neon photodisintegration 20 Ne(7, a) 16 0, and oxygen 
burning 16 0( 16 0,7) 32 S(7, a) 28 Si, all release a particles that can be captured by 16 0, 20 Ne, 
or 24 Mg through the reactions listed above. 

In accordance with the foregoing remarks, we need only follow changes in eight active 
isotopes, namely 1 H, 4 He, 12 C, 14 N, 16 0, 20 Ne, 24 Mg, and 28 Si. Thus, in eqs. ©-©, or 
in eqs. ( !9el) -(l9fl). the index j runs over the active isotopes, from 1 to 8. [The number of 
active isotopes may be changed, provided that the nuclear reaction network is modified 
accordingly.] Other isotopes, such as 40 Ca or 56 Fe, are regarded as inert: they contribute to 
the EOS, but their abundances do not change; in particular, they do not undergo convective 
mixing. For consistency, then, their abundances should be uniform throughout the initial 
stellar configuration, and so they will remain. 



3.4 Neutrino Losses 



Neutrino losses are according to 



Itohetal 



(119961), accounting for neutrino formation pro- 



cesses of pair annihilation, photo annihilation, plasma decay, bremsstrahlung and (option 
ally) r ecombination. There is also an option for using the older fitting formulae of 
(119671 ). 



Beaudet et al 



3.5 Convection (diffusive-convective mixing) 

So long as the radiative 'temperature gradient' dlnT/dliap , defined by 

kL p 



R 



does not exceed the adiabatic 



AncGm 4pn ' 

d\nT(p,s.X) 
dlnp 



(44) 



(45) 
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the actual gradient V , which appears in eq. (19c)) . is equal to Vr , and the diffusion coefficients 
o~j in eq. ( pel) are all zero: there is neither convective heat transport, nor any convective 
mixing. We hope that the specific entropy s in ( T4"5l) will not be confused with s = r 2 . 
Convection is taken to set in whenever 

AV = V R - V A > 0. (46) 

he actual gr adient V is then calculated in accordance with the mixing length recipe 



( IMihalaslll978l ): consider the non-dimensional (inverse) convective efficiency parameter 

l6V2a R T* 



(47) 



(£/H)^Qpc P Tv r e 

where <7r is Stefan's constant; I is the mixing length, which we take to be a constant multiple 

(of order unity) of the pressure scale height H = v^/g , where v 2 = p/p is the squared thermal 

speed and g = Gm/r 2 is the local acceleration of gravity; Q = — din p(p, T, X)/dlnT ; cp is 

the specific heat at constant pressure; and r e = np£ . Let x be the root of the cubic equation 

3 



— x A + x 2 + 2b'x = 1, 
4o 



where 



(48) 



(49) 



b' = 6/VAV, 
Then the actual gradient is given by 

V = V A + (£ 2 + 26V)AV. (50) 
It is readily seen that V — > Va as b' —>■ , and V — > Vr as b' —>■ oo . 

Convective mixing is taken to be due to diffusion in a gas of particles — representing the 



convective elements — moving at the convective speed 

■,1/2 



H 



XVq 



(51) 



( iMihalad Il978l ) , with the mean free path t . In such a gas the diffusion coefficient is ~ v c £ . 

But in eqs. (|9el)-( l9fl) the derivatives are with respect to mass, not radius. We therefore set 

the convective diffusion coefficients equal to 

/dm\ 2 



0-; 



V dr 7 



v c £ = (Artr 2 p) 2 v c e, (52) 

the same for all species j . 

The code sometimes runs into difficulties with the foregoing convective diffusion coeffi- 
cients. We therefore retain an option whereby the last formula is replaced by a much simpler 
one: 

'AV\ 2 



k, 



(53) 
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where k c < 1 is a numerical coefficient. Its purpose is to ensure that convective mixing does 
not occur too suddenly. The value of k c is related to the evolutionary time scale, and ranges 
from ~ 0.1 for low-mass stars to ~ 0.001 for massive ones. 

Finally, the code has an option for introducing convective overshoot. This is done in a 
rather schematic way: at each iteration, after determining the convective zones in accordance 
with the inequality AV > , we repeat the determination of the zone boundaries, this 
time with AV = Vr — Va + Vos > , where Vos is a small, positive constant. The 
temperature gradient V and the convective diffusion coefficient are then determined by the 
foregoing formulae, but with the ne w, augmented, A V . We do not attempt to fix Vos by 



any dependence on local conditions (IPols et al 



1995|). 



3.6 Mass Loss 



The stellar mass may change with time at a prescribed rate M, according to boundary con- 
dition ( fl8l) . This rate is generally taken to be a function of the stellar parameters M*, L+, R* . 
Over the years, several formulae have been suggested in the literature, each fitting observa- 
tions of stars in a particular evolutionary phase. We mention them briefly below, with the 
mass loss rate (MLR) in units of M & yr _1 . 



1. The earliest such expression is Reimers's formula (lReimerall975h . derived from obser- 
vations of RGB stars, 



Mi 



Reim 



-4 X lO^Ve 



1 AL 



(54) 



where the coefficient r]R eim lies between 0.3 and 3.0. 

2. A fit for early ty pe O and B stars, with somewhat modified powers of AL, L+, R± , is 



given by 



Lamerd (119811 ): 

M Lam 



-10 



-4.83 i 



/L + \ 1-42/^0.61,^ 

V 30 



-0.99 



(55) 



V10 3 7 V30 7 

3. A modific a tion o f Reimers's MLR, allowing for a superwind on the AGB, is given by 
Baud fe HabinJ (jl983( ): 

M env n 



M BU = M Reim x 



M P , 



(56) 



where M em) ,o is the envelope mass at the base of the AGB. 



4. Another variation on Reimers's MLR, simil ar to M\, arn , is given by lNieuwenhuijzen fe de Jager 



(If990l ) (subsequent paper to 



de Jager et al. 



Mi 



NDJ 



19881 . where M was given as a function of T e ff, L): 
9.63 x 10- 15 Li- 42 i£ 81 Mf 16 . (57) 
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5. The strong increase in mass- loss rate during the AGB stage is rendered by the MLR 



formulae of 



Blocker! dl99fih. w 



shock-driven winds by 



Bowen 



r h ich ar e based on an investigation of long period variables and 
( 119881 1 Blocker's MLR formula is: 

M B i = 4.83 x \^M Rei Mz 2 AMS L l\ (58) 
and a variant, Mb2, has Mzams replaced by M*. 

6. Yet another modification of Rei mers's formula, intended f or cool winds that are no t 



driven by molecules or dust, is given by 



Schroder fc Cuntzl (l2005luSchr5der fc Cuntzl (120071 ): 



sc 



~Vsc- 



3.5 



1 + 



9e 



4300#, 



(59) 



\A000KJ 

with rjsc = 8(±1) x 10~ 14 . Here two new factors are included, taking into account the 
dependence of chromospheric height on surface gravity and the dependence of the mechanical 
energy flux on the effective temperature. 

In applying any of the MLR expressions, instead of turning it on suddenly, we multiply 
it by a Fermi weight function 

F ^ R ^ = 1 _|_ e (Rthre S h-R*)/(0mR thre3h ) ' ( 60 ) 

where Rthresh is an MLR threshold radius, which we typically choose between 1 and 50. 
Its precise value is not important, so long as the MLR is negligible for R = Rthresh- As 
R* increases, F(R ir ) varies smoothly near Rthresh from to 1, over a width of 0.05R t hresh- 
This prevents an on-off situation, which can ruin the convergence of the iteration process by 
which the difference equations of §2.21 are solved. 

The question remains, which formula to use? The code includes an algorithm that iden- 
tifies the evolutionary stage of the stellar model by testing various parameters (such as 
luminosity, radius, composition profiles) and their rates of change. Therefore, one may pass 
- in a smooth manner — from one formula to another. In this work, we used ( 154"|) for the 
RGB and (158]) for later stages. The parameter ^R e im was taken progressively higher with 
increasing initial mass. The effect of Rthresh and r^Reim on the results will be briefly discussed 
in section §4.41 



4 EVOLUTION SEQUENCES 

Using the evolution code described in the previous section, we performed calculations over 
a wide range of initial stellar masses and metallicities. In the following sections we address 
representative results, outcome of continuous calculations that yield complete evolutionary 
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tracks, starting from an initial pre-main-sequence state and ending with either a cooling 
white dwarf (for initial masses below 9 M@), or core collapse of a configuration resembling a 
supernova progenitor (for higher initial masses). We use the following acronyms: MS - main 
sequence; ZAMS - sero-age main sequence; pre-MS - pre-main-sequence; RGB - red giant 
branch; HeF - helium flash; HB - horizontal branch; AGB - asymptotic giant branch; TP - 
thermal pulse; WD - white dwarf; HRD - Hertzsprung-Russell diagram. Central properties 
are denoted by subscript c. 



4.1 Solar Model 



We started from a 'pre-MS' configuration of 1M Q , of uniform composit ion Y = 0.29 and 



0.018 — the latter with a heavy element distribution according to 



Grevesse fc Noels 



(119931 ) — and a radius of 2.7 R Q . With a mixing-length to scaleheight ratio a = I /Hp = 
2.5, this configuration reached the ZAMS after 0.05 Gyr. At an age of 4.60 Gyr — which 
includes the 0.05 Gyr from pre-MS to ZAMS — the model reached a radius of 1.006 R & , 
a luminosity of 1.009 L & , and central characteristics T c = 15.59 x 10 6 K, p c = 2.453 x 
10 17 dyncm -2 , p c = 157.9 gr cm -3 . We regard this as a good match to the pres ent sun, and 



the central characteristics i n agreement 



1995 



Turcotte et al. 



1998 



Morel et al 



length recipe uses the constants of 



with those obtained by other codes (e.g.. 



Reiter et al. 



2000). It should, perhaps, be noted that our mixing- 
MihalaJ ( 1978 ). and our choice of a = 2.5 may correspond 



to different values for other choices of the constants. 

Fig. [3] shows the evolutionary track in the HRD, where the various phases are marked: 
from pre-MS, through MS, RGB and core HeF, settling into stable core He burning, con- 
tinuing through AGB and thermal pulses up to the last He shell flashes — where a strong 
flash occurs, followed by a weaker one — and ending with a cooling 0.55 M CO-WD. The 
durations of the MS, RGB and HB stages are 10 Gyr, 1.5 Gyr and 78 Myr, respectively. The 
maximum radius and luminosity — attained on the AGB after some 11.7 Gyr of evolution 
(from ZAMS) — are 1.46 x 10 2 R & and 2.81 x 10 3 L Q , respectively. The maximum temper- 
ature throughout the evolution, 2.09 x 10 s K, is attained off-center, at the tip of the AGB. 
We terminated the calculation with a final CO-WD of radius Rwd = 2.13 x 10~ 2 R & , a 
central pressure p c ,wD = 6.94 x 10 22 dyncm~ 2 , a central density p c ,wD = 1-84 x 10 6 grcm~ 3 
and a core temperature of ~ 75 million K. 
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4.2 The Effect of Metallicity 

The effect of metallicity on stellar evolution is illustrated by a series of calculations for a 
model of solar mass and (Z,Y) values of (0.0001,0.24), (0.001,0.24), (0.018,0.29), (0.05,0.30) 
and (0.1,0.30), other physical and numerical parameters remaining fixed. The results are 
presented in Fig. H] by complete, continuous tracks in the H-R diagram. We note that an 
increase in metallicity has a similar effect to a decrease in the initial stellar mass: luminosities 
are lower and the durations of evolutionary phases are longer. For example, the MS phase 
lasts up to over 3 times longer, when Z increases from 10 -4 to 0.1. This result is mostly 
the consequence of the dependence of opacity on composition; at a lower metallicity, the 
opacity decreases, the star is able to radiate away its energy with greater efficiency, the 
stellar luminosity is therefore higher and timescales are correspondingly shorter. 

Apart from the apparent shift of the evolutionary tracks in the H-R diagram, and the 
different timescales, metallicity also affects the final masses. For Mi = 1 M , a final mass of 
0.57 M & was obtained for the lowest metallicity {Z = 0.0001), and 0.52 M Q for the highest 
one [Z = 0.1), as compared with 0.55 M Q , obtained for solar metallicity — an overall spread 
of almost 10%. 

4.3 Canonical Evolution Sequences 

We consider Population I (Pop. I) and Population II (Pop. II) stars, adopting metallicities 
of Z = 0.01 and Z = 0.001, respectively, and initial masses in the range 0.25 — 9M , 
leading to cooling WDs. The complete evolutionary tracks are shown in the two panels of 
Fig. [51 Timescales and the final WD masses and composition are given in the accompanying 
Table 1. It should be noted that the MS and RGB durations as shown in the table depend 
on the definition of the MS-turnoff point and beginning of the RGB, which involves some 
arbitrariness. The criterion we use for the MS turnoff is as follows: let t\ be the time when X c 
has decreased below 10~ 6 ; let X\ = \ogT e ff(ti) and yi = logL(tx). The turnoff time t 2 is the 
earliest time for which the distance between the points [x 2 = \ogT e ff(t 2 ) , y 2 = log L(t 2 )} and 
[xi , yi] in the [log T e f f , log L] plane exceeds 0.1. Similar criteria are used for other transitions 
between evolutionary stages. Time scales depend strongly on composition, especially on Z, 
decreasing with decreasing Z. Given differences in composition adopted in different studies, 
as well as differences in criteria defining evolutionary stages, a precise comparison between 
models is difficult to achieve. Nevertheless, we find excellent agreement, for example, between 
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our low-mass Pop. II models and corresponding ones calculated by others: for the 0.8 M 
and 1 M models, we find tms — 1-44 x 10 10 yr and 6.02 x 10 9 y r, resp ectively (see Table 1), 



while for the same ma sses and metallicit y, 
and 6.06 x 10 9 yr, and 



Charbonnel et al. 



Charbonnel et al 



199^), after modifying the input physics, find 



( 119961 ) find 1.51 x 10 10 yr 



1.43 x 10 10 yr and 6.85 x 10 9 yr. For Pop. I mod els, the spread in initial Z is larger, yet 
our results are compatible with those obtained by ISerenelli fc Fukugital (120071 ) for a grid of 
stelar models with Z = 0.019. 

The evolutionary tracks end with a cooling WD. A He- WD is obtained for the lower 
initial masses, 0.25 < Mj < 0.50 M & . The transition to a CO- WD occurs between 0.50 and 
0.80 M Q , and the heavier ONeMg-dominated WDs are obtained for initial masses higher 
than ~ 8 M (the transition mass being higher for the Pop. I stars). It should be noted, 
however, that especially for the Pop. I stars, the transition mass for obtaining a CO- WD 
rather than a He- WD is strongly dependent on the mass-loss rate assumed. For example, 
for an initial mass of 0.80 M , slightly increasing the mass-loss rate may result either in 
a He- WD, when the threshold for core helium burning is not reached, or, in an Extreme 
Horizontal Branch (EHB) star, when a 'delayed' core HeF takes place. The production of 
such hot (blue) HB stars for relatively low initial masses (from around 0.80 to slightly over 
1 Mq) and for a range of metallicities will be addressed in a subsequent paper. 

For both populations, a violent ignition of helium takes place in the core (but usually 
off-center, because of neutrino cooling) at the tip of the first giant branch for masses in 
the range 0.80 — 2 M . This is the well-known core HeF. The transition between HeF and 
quiet He ignition occurs at an initial mass between 2 and 3 M , depending mainly on 
composition and mass-loss rate. During the flash, the peak nuclear energy generation rate 
is in the range 5 x 10 7 < L nuc ^ max < 5 x 10 9 L G , decreasing with increasing initial mass, 
due to a corresponding decrease in the degree of electron degeneracy of the core material. 
It is worth noting that the luminosity of the star during the flash is unaffected by what is 
taking place in the core, despite the huge nuclear luminosity, which surpasses the luminosity 
obtained at any evolutionary stage. The overall duration of the flash (when L nuc is in excess 
of, say, 10 5 L Q ) is of the order of a few years. We note that during this stage time steps are 
automatically reduced down to days, then hours and minutes. Once the flash is over, it will 
take some extra 10 3 to 10 5 years before the star settles into stable core He burning, the HB 
phase. 

The well-known thermal pulses that arise as a result of the double shell-burning instabil- 
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ity, are clearly seen in the evolutionary sequences during the final stages of the AGB. Fig. [6] 
shows a typical example for a Pop. II, 2 M® model. The thermal pulses in this example span 
about 7 x 10 5 yr, and clear trends are evident, such as the monotonic decrease in effective 
temperatures with advancing pulses, along with an increase in the radial extension of the 
photosphere, which reflect the asymptotic evolution towards the redder tip of the AGB. 
Also evident is the fact that the bulk of mass-loss takes place precisely during this short 
phase, with the mass dropping from 1.90 to 0.63 M & - almost its final value. The mass of 
the H-depleted core increases during this phase from 0.58 to 0.62 M ; the mass of the inner 
He-depleted core increases from 0.48 to 0.54 M . Since the He profile is not as steep as 
the H profile, the mass of the He-depleted core is a matter of definition: here 'He-depleted' 
means Y < 10~ 6 . Taking the core boundary at the mid-point of the He profile yields a final 
He-depleted core mass of O.6OM . We should note that the total number of pulses in each 
evolutionary sequence is largely determined by the mass-loss law adopted. 



4.4 Mass-Loss Laws and Initial-Final Mass Relationship (IFMR) 



Using the complete evolutionary tracks for the mass range of 0.8 to 9 M , for both popula- 
tions Z = 0.01 and Z = 0.001, as listed in Table 1, we obtain a theoretical IFMR, displayed 
in Fig. [7] (solid and dashed black lines). We increased the number of points by adding results 
for masses of 2.5 and 3.5 M , and for Z = 0.02 and Z = 0.005 and masses of 1,3 and 



5 M & (m arked in Fig . [7] by different symbols). Sim i 



puted bv lMeng et all (120071) an d bv 



obtained by 



Catalan et al- 



ar re lationships have been recently com- 



(120081 ). the latter including earlier results 



Dominguez et al.l (11999). A different and indep endent source for such a rela- 



tionship is provided by observations (e.g., 



Weidemann 



2000), mainly of st a r clus ters, which 



Ferrario et al. 



(120051 ) (based on 



lead to empirical or semi-empirical linear relations, such as 
open-cluster data for the range 2.5 — 6.5 M ) and others that will be mentioned below. 
The curves obtained here show that the IFMR may be divided into three regions with 



different s 
results of 



opes: 1. A modera te slope for Mj < 3 M , which coincides with the tabulated 
Weidemann! (120001 ) plotted in Fig. [7J 2. A steeper slope for 3 < Af< < 4 M . 



3. Again, a more gra dual incr e ase u ntil the top end. We note that the 'new relation' as 
displayed in Fig. 2 of iHerwigi (119951 ). meant to fit only the best determined stars of the 
Hyades and Pleiades clusters, has a very similar shape to our curves, only shifted upwards 
from our Pop. I curve by about 0.05 M . 
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The dependence on metallicity is app arent from th e diver gence o 



in agreement with the conclu sions of e.g. 



curves as plotted in Fig. 5 of 



Meng et all (120071 ) or the 



Catalan et al 



Mi < 2 M & (in agreement with e.g. 



2008 . 



the two curves in 



Dominguez et al.l (119991 ) 



Catalan et al. 



'he e ffect of metallicity is negligible for 



(120081 )). but it increas es towards higher 



Meng et al. 



2007 



initial masses: the curves diverge by > 0.13 M at the top end Mj > 7 M Q . 
reach a difference of up to 0.4 M in the final masses derived from different metallicities, 
their study covering a broad metallicity range: Z in between 0.0001 and 0.1. They also notice 
a minimum of the IFMR for Z = 0.04. 

Various semi-empirical linear fits have been derived over the last decade. A few examples 

are: 



Ferrario et al. 



(120051 ) (based on open-cluster data for the range 2.5 — 6.5 M ; claiming 
that the IFMR can be modelled by a mean relationship about which there exists some 
intrinsic scatter, and that they 'cannot justify the use of any but a linear relationship to 



model the cluster data'): 

Mf = (0.10038 ± 0.00518)71^ + 0.43443 ± 0.01467 



(61) 



Dobbie et al. 



(120061 ) (a linear fit to some 27 WDs, members of clusters such as the Hyades, 
Praesepe, M35, NGC2516 and the Pleiades, over initial-mass range of 2.7 — 6 M ): 

Mf = (0.133 ± 0.015)M, + 0.289 ± 0.051 (62) 



Williams! (120071 ) (claiming that the IFMR is both linear and without any metallicity 
dependence) : 

M f = (0.132 ± 0.017)Mi + 0.33 ± 0.07 (63) 

Although the relations obtained, as shown in Fig. 0, are quite far from linear, the closest 
linear fit that we can suggest, without using any artificial anchoring, is 

M f = 0.08343 *Mi + 0.47321 (64) 

which falls slightly above the upper (Pop. II) curve around the lower initial masses (1.5 — 
2.5 Mq), and below the lower (Pop. I) cu rve for higher interm ediate masses, around 5 M . 



This fit is very similar to the linear fit of 



Ferrario et al 



(120051 ) (shown in Fig. [7j) , although 



the latter is limited to the range 2.5 to 6.5 M . 

Clearly, the relation obtained represents the set of parameters assumed, mostly those 
related to the mass-loss recipe. The value of r/Reim used here was linearly increased from 0.4 at 
0.8 M Q to 3.0 at 9 M . A preliminary comparison that we performed, keeping all parameters 
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fixed and changing only mass-loss laws, indeed showed some differences in the final WD 
masses, with a spread of less than 10%. More precisely, for our solar model parameters (see 
§4.3p . setting rjR ei = 0.6, Rthresh = 50, the derived final WD masses were all in the range 
0.53 — 0.57 M Q (or between 0.51 — 0.56 for slightly higher mass-loss rates obtained by using 
VRei — 1-0, Rthresh = 10). Performing the same comparison for 3 M (Z = 0.01), but using 
rj Rei = 2.0, we found final WD masses to be in the range 0.61 — 0.67 M . 



4.5 Massive Stars 

We now briefly consider Pop. I massive stars of initial masses in the range 16 — 64 M , 
typically, SN progenitors. Since nucleosynthesis calculations are limited in our code, we 
cannot follow the evolution all the way to the collapse of an iron core. However, we come 
quite close to it. These massive stars go through advanced nuclear burning stages, until 
a core composed of the end-product of our nuclear reactions network is obtained. Core 
masses range monotonically from 2.4M for the 64M initial mass and 1.7M for the 16M 
initial mass. The core is enveloped by layers of different composition, the outermost being 
predominantly helium. Envelope masses depend strongly on the mass loss law assumed. 

The core contracts, becoming degenerate and unstable, since its mass exceeds the Chan- 
drasekhar limit. As contraction accelerates, temperatures rise to a few 10 10 K, where electron- 
positron pairs are created, which enhances the instability, lowering the adiabatic exponent. 
Pair production replaces iron photodisintegration as the mechanism leading to core col- 
lapse. Density profiles throughout the stars are shown in Fig. [81 The code crashes when 
the collapse approaches free-fall, with the adiabatic exponent very close to 4/3 throughout 
the core. Since this point is somewhat arbitrary, the curves representing stars of different 
initial masses do not exhibit a perfectly regular (monotonic) behaviour; this is so metimes 



the case for evolutionary track s or characteristics of massive stars in the late stages (lArnett 



1996 



Umeda fc Nomoto 



20081 ). resulting from the complexity of the processes taking part 
in them, and the related parameters and thresholds. We do not claim that these calculations 
shed light on pre-supernova evolution; rather, we mention them here as an example of the 
robustness of the code, which is capable of dealing with complex processes under critical 
conditions without failing. 

Finally, adding the results obtained for lower masses of Pop. I, described in Section 14.31 
we show in Fig. [9] evolutionary tracks of the stellar central points in the (log T, log p) plane, 
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exhibiting the branching off between stars that end their lives as WDs, and stars that go 
through advanced nuclear burning stages, ending their lives in dynamic core collapse. 



4.6 Non-Canonical Evolution 

The term 'non-canonical' refers to stars of unusual internal structure and composition. Such 
configurations may result from stellar mergers, where the merging stars may be MS stars, 
giants, compact stars or any combination of different types. Stellar mergers are probably 
the progenitors of blue straggler stars (BSS), found to exist in environments of high stellar 
density, such as globular clusters or the cores of open clusters. 

As already mentioned, the main reason for developing the evolution code presented here 
was the need for an efficient and fast tool that could be integrated into the MODEST (Mod- 
elling DEnse STellar systems) collaboration, combining dynamical N-body calculations with 
hydrodynamics — the colliding or merging of stars — and stellar evolution, for the simulat- 
ing of dense stellar environments. Whereas for normal stars, it is possible to construct and 
tabulate pre-computed evolutionary tracks for the use of MODEST calculations, merger 
products, having completely unpredictable configurations, must be evolved in situ. 

A non-canonical initial model will be the product of a hydrodynamic merger calculation, 
usually by smoothed particle hydrodynamics (SPH) methods. The first step in adapting 
such a model to quasi-static stellar evolution calculations is to obtain a hydr ostatically re- 



l axed configuration. This is achieved by applying the quasi-dynamic method of lRakavy et al. 



( 119671 ). Instead of eqs. flH)-([2J), consider the equations 

1 d Air 



p dm 3 



■r 



3 



(65) 



— = _4 7rr 2 MMiZl _ 9m (66) 

dr dm r 2 ' 

where r(m, r) is regarded as a function of the mass coordinate m and the quasi-time r , and 

p(p, s, Y) is determined by the EOS. The quasi-time has no physical meaning: its purpose 

is provide asymptotically (i.e. for r — > oo) a hydrostatic solution. Equation (1661) is called 

quasi-dynamic because the correct dynamic equation would have d 2 r/dt 2 — with t the true 

time — on its left-hand side. 

Let the boundary conditions be r = at the center, and p = at the surface. For a 

given distribution of entropy s(m) , and of the number fractions, collectively denoted by 
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Y(m) , and an initial distribution of radii r(m, 0) , the foregoing equations are to be solved 
for r(m, r) (and p(m, r) , and p(m, r) ). 

Since the entropy s and the composition Y are not varied, the (quasi) motion is adiabatic: 
du = —pd(l/p) . Multiplying fl66l) by dr/dr and integrating over the mass of the star yields, 
after an integration by parts, 

where 

E = £' '(u-^pjdm (68) 

is the total energy, internal and gravitational. Equation ( 1671) shows that the energy decreases 
with quasi-time. If — for the given entropy and composition distributions — a minimum of E 
exists, the solution of eqs. ( 165]) - (166 1) must lead to it, and the resulting structure, of stationary 
energy, will be hydrostatic. If, on the other hand, a minimum of E does not exist, the 
configuration is dynamically unstable: E will then decrease indefinitely. 

Thus, the quasi-dynamic method either leads to a hydrostatic structure, or else detects 
dynamical instability. It can be applied to any initial density distribution, even a uniform 
one. With the EOS 

p (p 7 s , Y) = Kp 1+ ^ (69) 

it can be used to construct a polytrope (dynamically unstable when n ^ 3), which may 
serve as an initial 'fully convective' protostellar model of uniform entropy and composition. 
Of course, 'solution' of (165]) - (166]) entails the replacement of the differential equations by 



impl icit difference equations, which are then solved by an iterative process (jRakavy et al 



19670 . 



As preliminary examples, we evolved merger products for three pairs of Pop. II [Z = 
0.001) low-mass parent stars. The parent stars were evolved by our code from some pre- 
MS initial configuration, to an age when the more massive star of each pair was almost at 
terminal MS age (TAMS), the less massive star of the pair being, of course, at an earlier 
stage on the MS. A pair of 0.85 and 0.60 M Q parent stars was evolved for 11 Gyr; a pair 
of 1.00 and 0.60 M© for 6 Gyr; and finally, a pair of 1.40 and 0.60 M for 1.5 Gyr. To 
calculate structures of the merger products for t he above pairs o f paren t stars, we used the 
MMAS ('make me a star', version 1.6) package of[ 



Lombardi et al. 



(2002), which produces ID 



models that approximate results of detailed SPH calculations. We chose to perform head-on 
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collisions (zero periastron separation), so that effects of rotation were absent. Each resulting 
merger product was incorporated as is into our code, and upon obtaining a hydrostatically 
relaxed configuration by the 'quai-dynamic method' as explained above, calculation of the 
evolution was initiated. 

It might be worthwhile to note the difference between the wa y we treat the merger - 



product and the way t he non-canonical evolution is initiated by iGlebbeek et al.l (120081 ). 



Glebbeek Polsl ( 120081 ). As explained in these papers, what the authors did was to start 



from a ZAMS model of the correct mass, evolve it until the central Xh equalled that of 
the merger product and then evolve it further with a fictitious energy production until its 
entropy profile equalled that of the merger product. This was done in steps, during which 
the composition was gradually adjusted to that of the merger product. This process resulted 
in a hydrostatic configuration that had the given mass and correct entropy and composition 
profiles. In contrast, what we did was to make use of the merger product exactly as obtained 
by the collision calculation and subject it to the quasi-dynamic method. 

Table I4~B1 lists some details of the colliding stars and the resulting mergers: t co i is the time 
of collision (age to which the parent stars were evolved); M merger is the mass of the merger 
product (slightly less than the sum of parent star masses, because some mass was lost in 
the merger process); Y c is the central He mass-fraction, tms is the remaining MS lifetime of 
the merger-product, whereas T MSfi0unter is the MS duration of the canonical counterpart - a 
normal ('canonical') star of initial mass equal to that of the merger-product. The central He 
mass fraction generally depends on the stages to which the parent stars have been evolved 
- how close to TAMS was the more massive parent star, and correspondingly, how much 
hydrogen did the less massive star of the pair managed to burn during its limited MS 
evolution. It should be noted, for instance, that the MS duration of the 1.88 M & merger- 
product exceeds that of the lower-mass 1.48 M & merger-product; this is due to the greater 
amount of central hydrogen in the more massive merger-product. 

Fig [10] shows evolutionary tracks on HRD of the three merger products (solid lines) 
(0.85 + 0.60, 1.00 + 0.60, 1.40 + 0.60 - top to bottom), while dashed lines represent evolu- 
tionary tracks of the canonical counterparts. The non-canonical models, possessing excess 
thermal energy right after the merging process, all begin by gravitational contraction before 
settling on the MS, where they spend the time required for burning the remaining central 
hydrogen. It is only during the MS and early-RGB phases that the non-canonical track dif- 
fers from that of the canonical one. The non-canonical track is shifted slightly upwards (to 



A Stellar Evolution Code for Calculating Complete Tracks 29 



Table 1. Collisions of low-mass MS parent stars - characteristics of the parent stars and merger-products. Masses are in solar 
units; MS durations are in years. 



Mi 


M 2 


tool (Gyr) 


^merger 


Y c 


TMS 


I'M S .counter 


0.85 


0.60 


11.0 


1.34 


0.96 


3.75e8 


2.05e9 


1.00 


0.60 


6.0 


1.48 


0.98 


2.84e8 


1.43e9 


1.40 


0.60 


1.5 


1.88 


0.88 


2.91e8 


6.10e8 



higher luminosity); the shift is growing with increasing mass (as is clearly apparent in the 
blow-up panels on the right). Except for the insignificant differences in the shape of the last 
shell flash while traversing the HRD from the AGB tip to the cooling WD curve, the tracks 
almost exactly overlap from RGB onwards. 

Fig. [TT] shows, as an example, composition profiles of the 1.40 + 0.60 M & merger, with 
comparison to the 1.88 M Q canonical counterpart at the point when the latter's Y c equals 
that of the initial state of the merger product. In the top panel the H and He profiles of the 
merger-product are plotted together with those of the 1.40 and 0.60 parent stars. Central 
hydrogen is almost completely depleted for the more massive parent star, which is very close 
to its TAMS; the low-mass parent star, at early stages of its MS evolution, still has a large 
fraction of hydrogen. 



Already a decade ago 



Sills et al 



(119971 ) began investigating evolutionary scenarios of col- 
lisionally merged stars, with the aim of examining possible formation channels and properties 
of blue straggler stars in globular clusters. They present results of evolutionary calculations 
for seven head-on collisions. Among their results, we find for instance a MS duration of 
3.74 x 10 8 yr for their 0.80 + 0.60 M Q merger; although details of the collision, including 
abundances, might not be exactly comparable, this result seems to be in very good agreement 
with our derived MS duration of 3.75 x 10 8 yr for our 0.85 + 0.60 M similar merger. 

As mentioned, more extensive evolut i onary calculations for collision products have re- 



cently been performed by 



Glebbeek et al. 



fbOQSluGlebbeek k Pols! (120081 ). Nowadays, several 



pro cedures for performin g calculations of stellar collisions, such as t he mentioned 



MMAS 



by 



Lombardi et al 



Gaburov et al. 



(I2008h 



(120021 ) or MMAMS ('make me a massive star') by 
are available. As illustrated by the foregoing three examples, our code is able to import and 
initiate evolution for merger-products created by either of the above procedures. In future, 
it will be interesting to study non-canonical evolution merger-products over a wider range 
of masses and initial compositions (outcomes of various combinations of the parent stars), 
as well as mergers involving other types of stars, such as compact objects — the merging of 
WD-MS or WD- WD. 
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SUMMARY 

We have developed a stellar evolution code that is capable of calculating full evolutionary- 
tracks without interruption or intervention. The implicit numerical scheme is based on si- 
multaneous solution of the thermodynamic and composition equations on an adaptive grid. 
Time steps are self-adjusting according to numerical as well as evolutionary time-scale cri- 
teria. The code was applied to a large variety of examples: full evolutionary tracks for stars 
of a wide range of masses and metallicities, and non-canonical stars obtained from stellar 
mergers. We believe that these examples of stellar evolution calculations demonstrate the 
efficiency and rubustness of our new code. We mention, in particular, the ability of the code 
to deal with the core He flash, thermal pulses, WD cooling, core collapse, as well as non- 
canonical configurations. We thus expect it to be useful in extensive parameter studies — of 
both stellar physics and initial properties of stellar models — as well as in simulations of 
stellar clusters. 
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Table 2. Timescales and final WD masses and compositions - Populations I and II canonical sequences 
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Figure 1. A Schematic representation of our set of 49 opacity tables — spanning a triangular shape in [X, Xqo\ space — in 
between which interpolations are performed for a given metallicity Z. The 7 open circles along the x-axis denote the 7 tables 
for zero CO excesses. Each point within the remaining 21 dots represents two tables: the excess being completely in carbon for 
one and completely in oxygen for the other (such as noted as example for the (X = .30 — Z/2, Xqq = -70 — Z/2) position). 
The hypotenuse of the triangle relates to zero helium mass fraction (Y = 0). 



This paper has been typeset from a TgX/ ETgX file prepared by the author. 
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Figure 2. Profiles of internal structure at three snapshots during evolution of a solar model - Mid-MS (left), tip of RGB 
(middle) and cooling WD (right). Top panels display profiles of opacity (solid black), density (dashed blue) and temperature 
(dot-dahsed red). Bottom panels display internal composition in terms of elemental mass fractions - hydrogen (solid red), 
helium (solid blue), and excesses of carbon and oxygen (dot-dashed magenta and cyan, respectively; values for the excesses 
are representing closely those of total C and O mass fractions). Dashed black plots display m/Mtot\ Mtot equal 1.00, 0.81 
and 0.55 Mq for the three profiles, left to right, respectively. Only in the rightmost panel are the carbon and oxygen excesses 
non-zero (post core helium burning). 




Figure 3. Solar model HRD - A complete evolutionary track as obtained for IMq, Y = 0.29, Z = 0.018 and mixing-length 
parameter a = 2.5. 
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Figure 4. Complete tracks on HRD for various metallicities Z = 0.0001 to 0.1 for IMq. MS effective temperatures and 
luminosities decrease with increasing metallicity; consequently - durations of MS increase (by a factor of over 3 from the lowest 
value of Z to the highest). 
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Figure 5. Complete tracks on HRD — Pre-MS to cooling WD — for initial (ZAMS) masses in the range 0.25 to 9.0 Mq. Top: 
population I models (Z = 0.01, Y = 0.28). Bottom: population II models (Z = 0.001, Y = 0.24). 
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Figure 6. Thermal pulses during TP-AGB for the 2 M Q Pop.II (Z = 0.001, Y = 0.24) model. Top: complete track on HRD; 
plotted in thick red is the TP-AGB phase, for which the bottom panels are plotted. Bottom: Evolution of various characteristics 
during the thermal pulses phase. 
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Figure 7. IFMR - Final vs. initial masses as obtained for both our Pop. I (solid black) and Pop. II (dashed black) evolutionary 
sequences, for initial m asses in the range .8 to 9 Mq. The solid blue line is a linear fit to all values (Pop. I and II). We show for 
comparison the r evised |Weide mannl ||2000| ) semi-empirical relationship (Mi in the range 1 to 7 Mq), as well as the empirical 
linear relation bv lFerrario et al .1 l|2005l ~ (Mi in the range 2.5 to 6.5 Mq). See text for details. 




Figure 8. Density profiles in massive stars (Pop. I) with collapsed cores (legend shows initial masses). 
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Figure 9. Evolution of the central stellar density and temperature for Pop. I (Z = 0.018) models in the range 0.25 — 64 Mq. 
Dotted line has a slope of 3 (as obtained for the log p c — log T c relation of hydrostatic equilibrium under ideal gas law). Nuclear 
burning phases are marked along the tracks. 
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Figure 10. Non-canonical evolution - evolutionary tracks on HRD of the three merger products (solid) of low-mass MS 
parent stars, with comparison to their canonical counterparts - 'normal' initial configurations of equal mass (dashed). Triangles 
denote starting points for both the canonical and non-canonical evolutionary tracks. Panels on the right are a blow-up of the 
(MS, early-RGB) regions; thicker solid and dashed sections denote the extent of MS evolution phase. 
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Figure 11. Left: Composition profiles (top: H,He; bottom: C,N,0) of the 1.40 + 0.60 Mq merger-product, as obtained from 
the Make Me A Star ver 1.6 package for a head-on collision (see text), right after the configuration has been hydrostatically 
relaxed by our code — ready to be evolved. Shown in thin lines at the top panel are the H and He profiles of the 1.40 Mq (red) 
and 0.60 Mq (blue) parent stars, evolved to an age of 1.5 Gyr. Right: similar profiles for the canonical counterpart - a 'normal' 
1.88 Mq star - when its central He mass fraction equals that of the merger product. 



